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ABSTRACT 

We present predictions for the evolution of the galaxy luminosity function, number counts and 
redshift distributions in the IR based on the ACDM cosmological model. We use the combined 
GALFORM semi-analytical galaxy formation model and GRASIL spectrophotometric code to 
compute galaxy SE Ds including the rep rocessing of radiation by dust. The model, which is 
the same as that in Bau gh et al. I d2005l) . assumes two different IMFs: a normal solar neigh- 
bourhood IMF for quiescent star formation in disks, and a very top-heavy IMF in starbursts 
triggered by galaxy mergers. We have shown previously that the top-heavy IMF seems to be 
necessary to explain the number counts of faint sub-mm galaxies. We compare the model 
with observational data from the Spitzer Space Telescope, with the model parameters fixed 
at values chosen before Spitzer data became available. We find that the model matches the 
observed evolution in the IR remarkably well over the whole range of wavelengths probed by 
Spitzer. In particular, the Spitzer data show that there is strong evolution in the mid-IR galaxy 
luminosity function over the redshift range z ~ — 2, and this is reproduced by our model 
without requiring any adjustment of parameters. On the other hand, a model with a normal 
IMF in starbursts predicts far too little evolution in the mid-IR luminosity function, and is 
therefore excluded. 

Key words: galaxies: evolution - galaxies: formation - galaxies: high-redshift - infrared: 
galaxies - ISM: dust, extinction 



1 INTRODUCTION 

In recent years, the evolution of galaxies at mid- and far-infrared 
wavelengths has been opened up for direct observational study by 
infrared telescopes in space. Already in the 1980s, the IRAS satel- 
lite surveyed the local universe in the IR, showing that much of 
present-day star formation is optically obscured, revealing a pop- 
ulation of luminous and ultra-luminous infrared galaxies (LIRGs 
with total IR luminosities L IR ~ 10 11 - 10 12 L Q and ULIRGs 
with Lir > 10 12 Lq), and providing the first hints of strong evo- 
lutio n in the number density of ULIRGs at recent cosmic epochs 
(e.g. IWright et al. II 1984 ISoifere? al. Ill987al : ISanders & Mirabel I 
1996). The next major advance came with the discovery by COBE 
of the cosmic far-IR background which has an en ergy density 
comparable to that in the optical/near-IR background ( Puge t"eF al. I 
ll996l ; |HauserefflZ.ll998l) . This implies that, over the history of the 
universe, as much energy has been emitted by dust in galaxies as 
reaches us directly in starlight, after dust extinction is taken into 
account. This discovery made apparent the need to understand the 
IR as much as the optical emission from galaxies in order to have 
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a complete picture of galaxy evolution. In particular, it is essen- 
tial to understand IR emission from dust in order to understand the 
cosmic history of star formation, since most of the radiation from 
young stars must have been absorbed by dust over the history of 
the universe, in ord er to account for the far-IR background (e.g. 
lHauseref a/ll998l) . 



Following these early discoveries, the ISO satellite enabled the 
first deep surveys of galaxies in the mid- and far-IR. The deep- 
est of these surveys were in the mid-IR at 15/rni, and probed the 
evolution of LIRGs and ULIRGs out to z ~ 1, showing strong 
evolution in these populations, and directly res olving most of the 
cosm i c infrared background a t that wavelength felbaz et al. II 19991 
120021 ; iGruppioni et al. 112002b. Deep ISO surv eys in the far-IR at 



jruppii 

170mn ( [Dole et al. I l200ll : IPatris et al. 1 120031) probed lower red- 
shifts, z ~ 0.5. Around the same time, sub-mm observations using 
the SCUBA instr ument on the JCMT revealed a huge population of 
high-z ULIRGs dSmail. Ivison & Blainll9l^ ; |Hugnes et al. 1 1998b 
which were s ubsequently found to h ave a redshift distribution peak- 
ing at z ~ 2 dChapman et al. ||2005|) , confirming the dramatic evo- 
lution in number density for this population seen at shorter wave- 
lengths and lower redshifts. The sub-mm galaxies have been stud- 
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ied in more detail i n subsequent SCUBA surveys (e.g. SHADES, 
iMortier et 0/120051) . 

Now observations using the Spitzer satellite dWerner et al. I 

|2004) . with its hugely increased sensitivity and mapping speed 
are revolutionizing our knowledge of galaxy evolution at IR wave- 
lengths from 3.6 to 160 pm. Spitzer surveys have allowed direct 
determinations of the evolution of the galaxy luminosity func- 
tion out to z ~ 1 in the rest-frame near-IR and to z ^ 2 in 
the mid-IR dLe Floc'h et al. 1 120051: IPerez-Gonzalez et al\ 120051 : 



B abb edge et al. 120061 : iFranceschini et al. I2006IK Individual galax - 
ies have been detected by Spitzer out to z ~ 6 dEvles et al. 112005b 
In the near future, the Herschel satellite (Pilbratt 2003) should 
make it possible to measure the far-IR luminosity function out to 
z ~ 2, and thus directly measure the total IR luminosities of galax- 
ies over most of the history of the universe. 

Accompanying these observational advances, various types of 
theoretical models have been developed to interpret or explain the 
observational data on galaxy evolution in the IR. We can distinguish 
three main classes of model: 

(a) Purely phenomenological models: In these models, 
the galaxy luminosity function and its evolution are described 
by a purely empirical expression, and this is combined with 
observationally-based templates for the IR spectral energy dis- 
tribution (SED). The free parameters in the expression for the 
luminosity function are then chosen to obtain the best match 
to some set of observational data, such as number counts 
and redshift distributions in different IR bands. These pa- 
rameters are purely descriptive and provide little insight into 
the physical processes wh ich control galaxy evolution . Exam - 
ples of these models are IPearson & Rowan-Robinson I ( 199€ ); 
Xu et al. I (fl998l); iBlain et al. I dl999l); IFranceschini et al.\ OOP lb; 
Charv & Elbaz (2001); Rowan-Robinson (2001); Lagache et al. 



d2003l) ; lGruppioni et al. I J2005h . 

(b) Hierarchical galaxy formation models with phenomeno- 
logical SEDs: In these models, the evolution of the luminosity func- 
tions of the stellar and total dust emission are calculated from a 
detailed model of galaxy formation based on the cold dark mat- 
ter (CDM) model of structure formation, including physical mod- 
elling of processes such as gas cooling and galaxy mergers. The 
stellar luminosity of a model galaxy is computed from its star 
formation history, and the stellar luminosity absorbed by dust, 
which equals the total IR luminosity emitted by dust, is calculated 
from this based on some treatment of dust extinction. However, 
the SED shapes required to calculate the distribution of the dust 
emission over wavelength from the total IR dust emission are ei- 
ther observationally-based temp lates (e.g.lGuiderdoni ~et al. II 199§l ; 
lDevriendT& Guiderd onTI l200d) or are purely phenomenological, 
e.g. a modified Planck function with an empirically chosen dust 
temperature (e.g. Kaviani et al. 2003). In this approach, the shape 
of the IR SED assumed for a model galaxy may be incompatible 
with its other predicted properties, such as its dust mass and radius. 

(c) Hierarchical galaxy formation models with theoretical 
SEDs: These models are similar to those of type (b), in that the 
evolution of the galaxy population is calculated from a detailed 
physical model of galaxy formation based on CDM, but instead 
of using phenomenological SEDs for the dust emission, the com- 
plete SED of each model galaxy, from the far-UV to the radio, is 
calculated by combining a theoretical stellar population synthesis 
model for the stellar emission with a theoretical radiative transfer 
and dust heating model to predict both the extinction of starlight 
by dust and the IR/sub-mm SED of the dust emission. The ad- 
vantages of this type of model are that it is completely ab initio, 



with the maximum possible theoretical self-consistency, and all of 
the model parameters relate directly to physical processes. For ex- 
ample, the typical dust temperature and the shape of the SED of 
dust emission depend on the stellar luminosity and the dust mass, 
and evolution in all of these quantities is computed self-consistently 
in this type of model. Following this modelling approach thus al- 
lows more rigorous testing of the predictions of physical models 
for galaxy formation against observational data at IR wavelengths, 
as well as shrinking the p arameter space of th e pr edictions. Ex- 
ample s of such models are lGranato et al. I d2000l) and lBaugh et al. I 
(2005). (An alternative modelling approach also based on theoret- 
ical IR SEDs but with a simplified tre atment of the assem bly of 
galaxies and halos has been presented by Granat oef al. I d2004l) and 
ISilvaefa/.ll2005h .) 

In this paper, we follow the third approach, with physical mod- 
elling both of galaxy formation and of the galaxy SEDs, includ- 
ing the effects of dust. This paper is the third in a series, where 
we c ombine the GALF ORM semi-analytical model of galaxy forma- 
tion dColeefa/.l200(jl) with the GRAS I L model for stellar and dust 
emission from galaxies dSilvaefaZ. Ill 998l) . The GALFORM model 
computes the evolution of galaxies in the framework of the ACDM 
model for structure formation, based on physical treatments of the 
assembly of dark matter halos by merging, gas cooling in halos, 
star formation and supernova feedback, galaxy mergers, and chem- 
ical enrichment. The GRAS I L model computes the SED of a model 
galaxy from the far-UV to the radio, based on theoretical models of 
stellar evolution and stellar atmospheres, radiative transfer through 
a two-phase dust medium to calculate both the dust extinction and 
dust heating, and a distribution of dust temperatures in each galaxy 
calculated from a detailed dus t grain model. In the first paper in 
the series dGranato et al. ||2000T) . we modelled the IR properties of 
galaxies in the local universe. While this model was very success- 
ful in explaining observations of the local universe, we found sub- 
sequently that it failed when confronted with observations of star- 
forming galaxies at high redshifts, predicting far too few sub-mm 
galaxies (SMGs) at z ~ 2 and Lyman-br eak galaxies (LBG s) at 
z ~ 3. Therefore, in the second paper (Baug h~ef al. Il2005h . we 
proposed a new version of the model which assumes a top-heavy 
IMF in starbursts (with slope x = 0, compared to Salpeter slope 
x — 1.35), but a normal solar neighbourhood IMF for quiescent 
star formation. In this new model, the star formation parameters 
were also changed to force more star formation to happen in bursts. 
This revised model agreed well with both the number counts and 
redshift distributions of SMGs detected at 850/im, and with the 
rest-frame far-UV luminosity function of LBGs at z ~ 3, while 
still maintaining consistency with galaxy properties in the local uni- 
verse such as the optical, near-IR and far-IR luminosity functions, 
and gas fractions, metallicities, morphologies and sizes. 

This same model of iBaugh et al. ] d2005h was found by 

iLe Delliou et and2005all2006^ to provide a good match to the ob- 
served evolution of the population of Lya-emitting galaxies over 
the redshift range z ~ 3 — 6. Support for the controversial assump- 
tion of a top-heavy IMF in burst s came from the s t udies o f chem- 
ical enrichment in this model by Nagashim a et al. I d2005al lbl). who 
found that the metallicities of both the intergalactic gas in galaxy 
clusters and the stars in elliptical galaxies were predicted to be sig- 
nificantly lower than observed values if a normal IMF was assumed 
for all star formation, but a greed much b etter if a top-heavy IMF in 
bursts was assumed, as in Baugh et al. . In this third paper in the 
series, we extend the Bau gh et al. I ( 120051) model to make predic- 
tions for galaxy evolution in the IR, and compare these predictions 
with observational data from Spitzer. We emphasize that all of the 
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model p arameters for the predictions presented in this paper were 
fixed by iBaugh et al. I prior to the publication of any results from 
Spitzer, and we have not tried to obtain a better fit to any of the 
Spitzer data by adjusting these parameter^. 

Our goals in this paper are to test our model of galaxy evo- 
lution with a top-heavy IMF in starbursts against observations of 
dust-obscured star-forming galaxies over the redshift range z ~ 
— 2, and also to test our predictions for the evolution of the stellar 
populations of galaxies against observational data in the rest-frame 
near- and mid-IR. The plan of the paper is as follows: In Section[2] 
we give an overview of the GALFORM and GRASIL models, fo- 
cusing on how the predictions we present later on are calculated. 
In Section [3] we compare the galaxy number counts predicted by 
our model with observational data in all 7 Spitzer bands, from 3.6 
to 160 (im. In Section [4] we investigate galaxy evolution in the 
IR in more detail, by comparing model predictions directly with 
galaxy luminosity functions constructed from Spitzer data. In Sec- 
tion [5] we present the predictions of our model for the evolution 
of the galaxy stellar mass function and star formation rate distribu- 
tion, and investigate the insight our model offers on how well stel- 
lar masses and star formation rates can be estimated from Spitzer 
data. We present our conclusions in Section|6] In the Appendix, we 
present model predictions for galaxy redshift distributions in the 
different Spitzer bands, to assist in interpreting data from different 
surveys. 



2 MODEL 

In this paper use the GALFORM semi-analytical model to predict 
the physical properties of the galaxy population at different red- 
shifts, and combine it with the GRASIL spectrophotometric model 
to predict the detailed SEDs of model galaxies. Both GALFORM 
and GRASIL have been described in detail in various previous 
papers, but since the descriptions of the different model compo- 
nents, as well as of our particular choice of parameters, are spread 
among different papers, we give an overview of both of these here. 
GALFORM is described in fl2~T1 and GRASIL in Particularly 
important features of our model are the triggering of starbursts by 
mergers (discussed in §2. 1.4t and the assumption of a top-heavy 
IMF in starbursts (discussed in £]2.1.7t . We further discuss the 
choice of model parameters in £]2.3I Readers who are already fa- 
miliar with the Baug h~er al. I J2005t) model can skip straight to the 
results, starting in Sj3] 

2.1 GALFORM galaxy formation model 

We compute the formation and evolution of galaxies within the 
framework of the ACDM model of structure formation using the 
semi-analytical galaxy formation model GALFORM. The general 
methodology and approximations behind the GALFORM m odel are 
set out in detail in lCole et al. ll l200bT) (see also the review by Baugh 
(2006)). In summary, the GALFORM model follows the main pro- 
cesses which shape the formation and evolution of galaxies. These 
include: (i) the collapse and merging of dark matter halos; (ii) the 

1 A closely related model of galaxy formation obtained by applying 
GALFORM principles to the Mille nnium simulat ion of Springel et al. IfeOOa) 
has recently been published by iBower et al. ] 120031 This model differs 
from the current one primarily in that it includes feedback from AGN activ- 
ity, but does not have a top-heavy IMF in bursts. We plan to investigate the 
IR predictions of this alternative model in a subsequent paper. 



shock-heating and radiative cooling of gas inside dark halos, lead- 
ing to the formation of galaxy disks; (iii) quiescent star formation 
in galaxy disks; (iv) feedback both from supernova explosions and 
from photo-ionization of the IGM; (v) chemical enrichment of the 
stars and gas; (vi) galaxy mergers driven by dynamical friction 
within common dark matter halos, leading to the formation of stel- 
lar spheroids, and also triggering bursts of star formation. The end 
product of the calculations is a prediction of the numbers and prop- 
erties of galaxies that reside within dark matter haloes of different 
masses. The model predicts the stellar and cold gas masses of the 
galaxies, along with their star formation and merger histories, their 
sizes and metallicities. 

The prescriptions and parameters for the different processes 
which we use in t his paper are identical to those adopted by 
Bau gh et all d2005h . sut differ in several important respects from 
Cole et ai] ( 2000). All of these parameters were chosen by com- 
parison with pie-Spitzer observational data. The background cos- 
mology is a spatially flat CDM universe with a cosmological con- 
stant, with "concordance" parameters Q, m = 0.3, Q,a = 0.7, 
tt b = 0.04, and h = Ho/(100kms -1 Mpc _1 ) = 0.7. The am- 
plitude of the initial spectrum of density fluctuations is set by the 
r.m.s. linear fluctuation in a sphere of radius 8/i~ 1 Mpc, erg = 0.93. 
For completeness, we now summarize the prescriptions and param- 
eters used, but gi ve details mainly where they differ from those in 
ICole et al. I (2000), or where they are particularly relevant to pre- 
dicting IR emission from dust. 



2.1.1 Halo assembly histories 

As in lCole et al. we describe the assembly histories of dark 

matter halos through halo merger trees which are calculated using 
a Monte Car lo method based on t he extended Press-Schechter ap- 
proach (e.g. lLacev &Cole|[l993T) . The process of galaxy forma- 
tion is then calculated separately for each halo merger tree, follow- 
ing the baryonic physi cs in all of the sepa rate branches of the tree. 
As has been shown bv lHellve? aU hoOD, the statistical properties 
of galaxies calculated in semi-analytical models using these Monte 
Carlo merger trees are very similar to those computed using merger 
trees extracted directly from N-body simulations. 



2.1.2 Gas cooling in halos 

The cooling of gas i n halos is c alculated using the same simple 
spherical model as in lCole et al. (2000). The diffuse gas in halos 
(consisting of all of the gas which has not previously condensed 
into galaxies) is assumed to be shock-heated to the halo virial tem- 
perature when the halo is assembled, and then to cool radiatively by 
atomic processes. The cooling time depends on radius through the 
gas density profile, which is assumed to have a core radius which 
grows as gas is removed from the diffuse phase by condensing into 
galaxies. The gas at some radius r in the halo then cools and col- 
lapses to the halo centre on a timescale which is the larger of the 
cooling time t coo \ and the free-fall time tg at that radius. Thus, for 
tcooi(f) > tg(r), we have hot accretion, and for t ool(f) < ts(r), 
we have cold accretion^ In our model, gas only accretes onto the 
central galaxy in a halo, not onto any satellite galaxies which share 

2 Note that contrary to claims by iBirnboim & Dekel I Jlobll) . the process 
of "cold acc retion", if not t he nam e, has always been part of semi-analytical 
models (see Croton et al. ( 2006) for a detailed discussion) 
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that halo. We denote all of the diffuse gas in halos as "hot", and all 
of the gas which has condensed into galaxies as "cold". 

2.1.3 Star formation timescale in disks 

The global rate of star formation ip in galaxy disks is assumed to be 
related to the cold gas mass, M gas , by ij) = Mgas/r^disk, where 
the star formation timescale is taken to be 

T^disk^r^K^OOkms- 1 )"*, (1) 

where V c is the circular velocity of the galaxy disk (at its half- 
mass radius) and r*o is a constant. We adopt values r»o = 8 Gyr 
and a* = —3, chosen to reproduce the observed relation between 
gas mass an d B-band luminosity for present-day disk galaxies. As 
discussed in lBaugh et al. I J2005h . this assumption means that the 
disk star formation timescale is independent of redshift (at a given 
V c ), resulting in disks at high redshift that are much more gas-rich 
than at low redshift, and have more gas available for star formation 
in bursts triggered by galaxy mergers at high redshift. 

2.1.4 Galaxy mergers and triggering of starbursts 

In the model, all galaxies originate as central galaxies in some halo, 
but can then become satellite galaxies if their host halo merges into 
another halo. Mergers can then occur between satellite and central 
galaxies within the same halo, after dynamical friction has caused 
the satellite galaxy to sink to the centre of the halo. Galaxy mergers 
can produce changes in galaxy morphology and trigger bursts. We 
classify galaxy mergers according to the ratio of masses (including 
stars and gas) M2/M1 ^ 1 of the secondary to primary galaxy 
involved. We define mergers to be major or m inor according to 
wheth er M2/M1 > / c m p or M 2 /M-i < / cmp (Kauffman n et al. I 
1 1993b . In major mergers, any stellar disks in either the primary or 
secondary are assumed to be disrupted, and the stars rearranged into 
a spheroid. In minor mergers, the stellar disk in the primary galaxy 
is assumed to remain intact, while all of the stars in the secondary 
are assumed to be added to the spheroid of the primary. We adopt a 
threshold / C ni p = 0.3 for maj or mergers, co nsistent with the results 
of numerical simulations (e.g. lBarnes 11 998). which reproduces the 
observed present-day fraction of spheroidal galaxies. We assume 
that major mergers always trigger a starburst if any gas is present. 
We also assume that minor mergers can trigger bursts, if they sat- 
isfy both M2 /Mi > /burst and the gas frac tion in the disk of th e 
primary galaxy exceeds / gas ,crit. Following Bau gh et al. 1 d2005t) . 
we adopt /burst = 0.05 and / ga s,crit = 0.75. The parameters for 
bursts in minor mergers were motivated by trying to explain the 
number of sub-mm galaxies. An important consequence of assum- 
ing eqn.lQJ for the star formation timescale in disks, combined with 
the triggering of starbursts in minor mergers, is that the global star 
formation rate at high redshifts is dominated by bursts , while that at 
low redshifts it is dominated by quiescent disks (see Bau gh et al. I 
for a detailed discussion of these points). 

In either kind of starburst, we assume that the burst consumes 
all of the cold gas in the two galaxies involved in the merger, and 
that the stars produced are added to the spheroid of the merger rem- 
nant. During the burst, we assume that star formation proceeds ac- 
cording to the relation ip — M gaa /r*, burst - For the burst timescale, 
we assume 

T*, burst — iriciX [/dyn7"dyn,sph , 7~* .burst ,min] , (2) 

where rd yn ,sph is the dynamical time in the newly-formed spheroid. 



We adopt /d yn = 5 and r«, burst . mm = 0.2 Gyr (these parame- 
ters were chosen bv lBaugh et al. I J2005t) to allow a simultaneous 
match to the sub-mm number counts and to the local 60/im lu- 
minosity function). The star formation rate in a burst thus decays 
exponentially with time after the galaxy merger. It is assumed to be 
truncated after 3 e-folding times (where the e -folding time takes ac- 
count of stellar recycling and feedback - see lGranato et al. I fcOOOl) 
for details), with the remaining gas being ejected into the galaxy 
halo at that time. 

2. 1.5 Feedback from photo-ionization 

After the intergalactic medium (IGM) has been reionized at redshift 
Zrcion, the formation of low-mass galaxies is inhibited, both by the 
effect of the IGM pressure inhibiting collapse of gas into halos, and 
by the reduction of gas cooling in halos due to the photo-ionizing 
background. We model this in a simple way, by assuming that for 
z < Zrcion, cooling of gas is completely suppressed in halos with 
circular velocities V c < Krit - We adopt Vcr i t = 60 kins -1 , based 
on the detailed modelling bv lBenson et al. I d2002h . We assume in 
this paper that reionizatio n occurs at z Ic i on = 6, for consistency 
with iBaugh et al. I d2005h . but increasing this to z r cion ~ 10 in 
line with the WMAP 3 -year estimate of th e polarization of the mi- 
crowave background (Spergere7aT|[2006) has no significant effect 
on the model results presented in this paper. 

2.1.6 Feedback from supernovae 

Photo-ionization feedback only affects very low mass galaxies. 
More important for most galaxies is feedback from supernova ex- 
plosions. We assume that energy input from supernovae causes gas 
to be ejected from galaxies at a rate 

M ci = p(Vc) i> = [fteh(vy + /MK)] i> 0) 

The supernova feedback is assumed to operate for both quiescent 
star formation in disks and for starbursts triggered by galaxy merg- 
ers, with the only difference being that we take V c to be the circular 
velocity at the half mass radius of the disk in the former case, and 
at the half-mass radius of the spheroid in the latter case. For sim- 
plicity, we keep the same feedback parameters for starbursts as for 
quiescent star formation. 

The supernova feedback has two components: the reheating 
term /3 TC hip describes gas which is reheated and ejected into the 
galaxy halo, from where it is allowed to cool again after the halo 
mass has doubled through hierachical mass build-up. For this, we 
use the parametrization of ICole et 

/3reh = (K/Vhotr Qhot , (4) 

where we adopt parameter values Vhot = 300 kms -1 and a n ot = 
2. The reheating term has the largest effect on low-mass galaxies, 
for which ejection of gas from galaxies flattens the faint-end slope 
of the galaxy luminosity function. 

The second term /3 sv ,tp in eqn.lO is the superwind term, which 
describes ejection of gas out of the halo rather than just the galaxy. 
Once ejected, this gas is assumed never to re-accrete onto any halo. 
We model the superwind ejection efficiency as 

/3 S „ = f„ min [1, (K/Vs„r 2 ] (5) 

based on lBenson et al. Id200"3l). We adopt pa r amete r values / sw = 2 
and Vsw = 200 kms -1 , as in lBaugh et al. I J2005I) . The superwind 
term mainly affects higher mass galaxies, where the ejection of gas 
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from halos causes an increase in the cooling time of gas in halos 
by reducing the gas densities. This brings the predicted break at 
the bright end of the local galaxy lumi nosity function into agree- 
ment with observations, as discussed in lBenson et al. I d2003h ■ The 
various parameters for supernova feedback are thus chosen in or- 
der to match the observed present-day optical and near-IR galaxy 
luminosity functions, as well as the galaxy metallicity-luminosity 
relation. 

We note that the galaxy formation model in this paper, un- 
like some other recent semi-analytical models, does not include 
AGN feedback. Instead, the role of AGN feedback in reducing the 
amount of gas cooling to form massive galaxies is taken by su- 
perwinds driven by supernova explosions. T he first semi-analytica l 
model to include AGN feedback was that of iGranato et al. I d2004l) . 
who introduced a detailed model of feedback from QSO winds 
during the formation phase of supermassive black holes (SMBHs), 
with the aim of explaining the co-evolution of the spheroidal com- 
ponents of ga laxies and their SMBHs. The predictions of the 
IGranato et al. I model for number co unts and redshift d istributions 
in the IR have been computed by Isilva et oTI (2005) using the 
GRASIL spectrophotometri c model, and compa red to ISO and 
Spitzer data. However, the IGranato et al. ] d2004l) model has the 
limitations that it does not include the merging of galaxies or of 
dark halos, and does not treat the formation and evolution of galac- 
tic disks. More recently, several semi-analytical models have been 
published which propose that heating of halo gas by relativistic 
jets from an AGN in an optically inconspicuous or "radio" mode 
can balance radiative cooling of gas in high- mass halos, thus sup- 



pressing hot accret i on of gas onto galaxi e s dBower et al. 
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2006; 
2007b . 



However, these AGN feedback models differ in detail, and all are 
fairly schematic. None of these models has been shown to repro- 
duce the observed number counts and redshifts of the faint sub-mm 
galaxies. 

The effects of our superwind feedback are qualitatively quite 
similar to those of the radio-mode AGN feedback. Both superwind 
and AGN feedback models contain free parameters, which are ad- 
justed in order to make the model fit the bright end of the ob- 
served present-day galaxy luminosity function at optical and near- 
IR wavelengths. However, since the physical mechanisms are dif- 
ferent, they make different predictions for how the galaxy lumi- 
nosity function should evolve with redshift. Current models for the 
radio-mode AGN feedback are very schematic, but they have the 
advantage over the superwind model that the energetic constraints 
are greatly relaxed, since accretion onto black holes can convert 
mass into energy with a much higher efficiency than can supernova 
explosions. We will investigate the predictions of models with AGN 
feedback for the IR and sub-mm evolution of galaxies in a future 
paper. 

2.1.7 The Stellar Initial Mass Function and Chemical Evolution 

Stars in our model are assumed to form with different Initial Mass 
Functions (IMFs), depending on whether they form in disks or in 
bursts. Both IMFs are taken to be piecewise power laws, with slopes 
x defined by AN /A In m <x m~ x , with N the number of stars and 
m the stellar mass (so the Salpeter slope is x — 1.35), and covering 
a stellar mass range 0.15 < m < 12OM0. Quiescent star forma- 
tion in galaxy disks i s assumed to have a solar neighbourhood IMF, 
for which we use the lKennicutt 1 dl983h paramerization, with slope 
x = .4 for m < M and x = 1.5 for m > M e . (The lKennicutt I 
dl983l) IMF is similar to other popular parametrizations of the solar 



neighbourhood IMF, such as that o f lKroupald200ll) .) Bursts of star 
formation triggered by galaxy mergers are assumed to form stars 
w ith a top-heavy IMF with slope x — 0. As discussed in detail 
in lBaugh etal. 1 d2005l) . the top-heavy IMF in bursts was found to 
be required in order to reproduce the observed number counts and 
redshift distr ibutions of the faint sub-m m galaxies. Furthermore, 
as shown by iNagashima et al. 1 d2005al lbh. the predicted chemical 
abundances of the X-ray emitting gas in galaxy clusters and of the 
stars in elliptical galaxies also agree better with observational data 
in a model with the top-heavy IMF in bursts, rather than a universal 
solar neighbourhood IMF. 

A variety of other observational evidence has accumulated 
which suggests that the IMF in some environments may be top- 
heavy compared to the solar neighbourhood IMF. iRieke et al. I 
d 19931) argued for a top-heavy IMF in the nearby star burst M82, 
based on modelling its integrated properties, while Par ra et al\ 
(2007) found possible evidence for a top-heavy IMF in the ultra- 
luminous starburst Arp220 from the relative numbers of super- 
novae of different types observed at radio wavelengths. Evidence 
has been found for a top-heavy IMF in some st ar clusters in in- 
tense ly star-forming regions, both in M82 (e. g. McCradv et al. 
2003), and in our own Galaxy (e.g. lFiger et al. II 1 9991 ; Is tolte et al. 
2005; Haravam a et al. 1 120071) . Observations of both the old and 
young stellar populations i n the central 1 pc of our Galaxy 
also favour a top-heavy IM F dPaumard et al. luOOrj ; iManess et al. I 
l2007h . lFardal et all (2006) found that reconciling measurements of 
the optical and IR extragalactic background with measurements of 
the cosmic star formation history also seeme d to require an average 
IMF that was somewhat top-heavy. Finally, Ivan Dokkum (2007) 
found that reconciling the colour and luminosity evolution o f early- 
type g alaxies in clusters also favoured a top-heavy IMF. lLarson I 
( 1998) summarized other evidence for a top-heavy IMF during the 
earlier phases of galaxy evolution, and argued that this could be 
a natural consequence of the temperature-depen dence of the J eans 
mass for gravitational instability in gas clouds. lLarson I ( 120051) ex- 
tended this to argue that a top-heavy IMF might also be expected 
in starburst regions, where there is strong heating of the dust by the 
young stars. 

In our model, the fraction of sta r formation occurin g in the 
burst mode increases with redshift (see Bau gh etal. Id2005f) ), so the 
average IMF with which stars are being formed shifts from being 
close to a solar neighbourhood IMF at the present-day to being very 
top-heavy at high redshift. In this model, 30% of star formation 
occured in the burst mode when integrated over the past history of 
the universe, but only 7% of the current stellar mass was formed 
in bursts, because of the much larger fraction of mass recycled by 
dying stars for the top-heavy IMF. We note that our predictions for 
the IR and sub-mm luminosities of starbursts are not sensitive to 
the precise form of the top-heavy IMF, but simply require a larger 
fraction of m ~ 5 — 2QMq stars relative to a solar neighbourhood 
IMF. 

In this paper, we calculate chemical evolution using the instan- 
taneous recycling approximation, which depends on the total frac- 
tion of mass recycled from dying stars (R), and the total yield of 
heavy elements (p). Both of these parameters depend on the IMF. 
We use the results of stellar evolution computa tions to calculate 
values o f R and p consistent with each IMF (see Nag ashima et al. I 
(2005a) for details of the stellar evolution data used). Thus, we use 
R = 0.41 and p = 0.023 for the quiescent IMF, and R = 0.91 and 
p — 0.15 for the burst IMF. Our chemical evolution model then 
predicts the masses and total metallicities of the gas and stars in 
each galaxy as a function of time. 
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2.1.8 Galaxy sizes and dust masses 

For calculating the extinction and emission by dust, it is essential to 
have an accurate calculation of the dust optical depths in the model 
galaxies, which in turn depends on the mass of dust and the size 
of the galaxy. The dust mass is calculated from the gas mass and 
metallicity predicted by the chemical enrichment model, assuming 
that the dust-to-gas ratio is proportional to metallicity, normalized 
to match the local ISM value at solar metal l icity. T he sizes of galax- 
ies are computed exactly as in lCole et al. I d2000t> : gas which cools 
in a halo is assumed to conserve its angular momentum as it col- 
lapses, forming a rotationally-supported galaxy disk; the radius of 
this disk is then calculated from its angular momentum, includ- 
ing the gravity of the disk, spheroid (if any) and dark halo. Galaxy 
spheroids are built up both from pre-existing stars in galaxy merg- 
ers, and from the stars formed in bursts triggered by these mergers; 
the radii of spheroids formed in mergers are computed using an 
energy conservation argument. In calculating the sizes of disks and 
spheroids, we include the adiabatic contraction of the dark halo due 
to the gravi ty of the ba r yonic components. This m odel was tested 
for disks bv lCole et al. I d2000h and for spheroids bv lAlmeida et al. I 
(2007) (see also Coenda et al. in preparation, and Gonzalez et al. 
in preparation). During a burst, we assume that the gas and stars 
involved in the burst have a distribution with the s ame half-mass 
radius as the spheroid (i.e. r) — 1 in the notation of iGranato et al. I 
(2000), who used a value r\ = 0.1). 



2.2 GRAS I L model for stellar and dust emission 

For each galaxy in our model, we compute the spectral energy dis- 
tribut i on using the spectrop hotometric model GRAS I L dSilva et al] 
1 19981 ; IGranato et al. Il2000h . GRAS I L computes the emission from 
the stellar population, the absorption and emission of radiation by 
dust, and also ra dio emission (therma l and synchrotron) powered 
by massive stars l lBressan et al. I2002T) . 



with a continuous distribution of grain sizes (varying between 8A 
and 0.25 fim), and also Polycyclic Aromatic Hydrocarbon (PAH) 
molecules with a distribution of sizes. The equilibrium temperature 
in the local interstellar radiation field is calculated for each type 
and size of grain, at each point in the galaxy, and this information 
is then used to calculate the emission from each grain. In the case of 
very small grains and PAH molecules, temperature fluctuations are 
important, and the probability distribution of the temperature is cal- 
culated. The detailed spectrum o f the PAH emission i s obtained us- 
in g the PAH cross-se ctions from lLi & Draine I d200lh . as described 



Vega et al. U2005I) . The grain size distribution is chosen to match 



the mean dust extinction curve and emissivity in the local ISM, and 
is not varied, except that the PAH abundance in mole cular clouds 
is ass umed to be 10 -3 of that in the diffuse medium ( Veg a et al\ 
l2005h . 

(vi) Radio emission from ionized gas in HII regions and from syn- 
chrotron radiation from relativistic electrons accele rated in super- 
nova r emnant shocks are calculated as described in Bressa n et al. I 
j2002h . 

The output from GRAS I L is then the complete SED of a 
galaxy from the far-UV to the radio (wavelengths lOOA < A < 
lm). The SED of the dust emission is computed as a sum over the 
different types of grains, having different temperatures depending 
on their size and their position in the galaxy. The dust SED is thus 
intrinsically multi-temperature. GRAS I L has been shown to give an 
excellent match to the measured S EDs of both quiescent (e.g. M51) 
and starburst (e.g. M82) galaxies dSilva et al. 1 199 8; Bressa n et al. I 
l2002h . 

The assumption of axisymmetry in GRAS I L is a limitation 
when considering starbursts triggered by galaxy mergers. However, 
observations of local ULIRGs imply that most of the star formation 
happens in a single burst component after the galaxy merger is sub- 
stantially complete, so the assumption of axisymmetry for the burst 
component may not be so bad. 



2.2.1 SED model 

The main features of the GRAS I L model are as follows: 

(i) The stars are assumed to have an axisymmetric distribution in a 
disk and a bulge. Given the distribution of stars in age and metal- 
licity (obtained from the star formation and chemical enrichment 
history), the SED of the stellar population is calculated using a 
population synthesis model based on t he Padova stellar evo lution 
tracks and Kurucz model atmospheres JBressan et al. 1 11998). This 
is done separately for the disk and bulge. 

(ii) The cold gas and dust in a galaxy are assumed to be in a 2-phase 
medium, consisting of dense gas in giant molecular clouds embed- 
ded in a lower-density diffuse component. In a quiescent galaxy, 
the dust and gas are assumed to be confined to the disk, while for 
a galaxy undergoing a burst, the dust and gas are confined to the 
spheroidal burst component. 

(iii) Stars are assumed to be born inside molecular clouds, and then 
to leak out into the diffuse medium on a timescale t csc . As a result, 
the youngest and most massive stars are concentrated in the dustiest 
regions, so they experience larger dust extinctions than older, typ- 
ically lower-mass stars, and dust in the clouds is also much more 
strongly heated than dust in the diffuse medium. 

(iv) The extinction of the starlight by dust is computed using a ra- 
diative transfer code; this is used also to compute the intensity of 
the stellar radiation field heating the dust at each point in a galaxy. 

(v) The dust is modelled as a mixture of graphite and silicate grains 



2.2.2 GRASIL parameters 

The main parameters in the GRASIL dust model are the fraction 
/ mc of the cold gas which is in molecular clouds, the timescale t csc 
for newly-formed stars to escape from their parent molecular cloud, 
and the cloud masses M c and radii r c in the combination M c jr\, 
which determines the dust optical depth of the clou ds. We assume 
L = 0.25, M c = 10 6 M Q and r c = 16pc as in IGranato et al.\ 
(2000), and also adopt the same geometrical parameters as in that 
paper. We make the followin g two changes in GRASIL parame - 
ters relative to lGranato et al. L as discussed in iBaughef al. U2005I) : 
(a) W e assume t BSC = 1 Myr in both disks and bursts (instead 
of the IGranato et al. I values t CBC = 2 and 10 Myr respectively). 
This value was chosen in order to obtain a better match of the pre- 
dicted rest-frame far-UV luminosity function of galaxies at z ~ 3 
to that measured for Lyman-break galaxies, (b) The dust emissivity 
law in bursts at long wavelengths is modified from e v oc v~ 2 to 
e„ oc v~ 1,5 for A > 100/im. This was done in order to improve 
slightly the fit of the model to the observed sub-mm number counts. 
In applyin g GRASIL to mo del the SEDs of a sample of nearby 
galaxies, Isiiva et al. I < fl998l) found that a similar modification (to 
e„ oc u" 1 ' 6 ) seemed to be required in the case of Arp220 (the only 
ultra-luminous starburst in their sample), in order to reproduce the 
observed sub-mm data for that galaxy. This modification in fact has 
little effect on the IR predictions presented in the pre sent paper, but 
we retain it for consistency with lBaugh et al. I d2005h . 
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2.2.3 Interface with GALFORM 

For calculating the statistical properties of the galaxy population 
from the combined GALFOR M+GRASIL model, w e follow the 
same strategy as described in iGranato et al. I d2000h . We first run 
the GALFORM code to generate a large catalogue of model galaxies 
at any redshift, and then run the GRASIL code on subsamples of 
these. For the quiescent galaxies, we select a subsample which has 
equal numbers of galaxies in equal logarithmic bins of stellar mass, 
while for the bursting galaxies, we select a subsample with equal 
numbers of galaxies in equal logarithmic bins of burst mass. For the 
burst sample, we compute SEDs at several different representative 
stages in the burst evolution, while for the quiescent sample, we 
only compute SEDs at a single epoch. Using this sampling strat- 
egy, we obtain a good coverage of all the different masses, types 
and evolutionary stages of galaxies, while minimizing the compu- 
tational cost of running the GRAS I L code. The statistical properties 
of the galaxy population are then obtained by assigning the model 
galaxies appropriate weights depending on their predicted number 
density in a representative cosmological volume. 

The outputs from the GALFORM galaxy formation model re- 
quired by GRASIL to calculate the galaxy SEDs are: the combined 
star formation history and metallicity distribution for the disk and 
bulge, the radii of both components, and the total mass of dust. The 
dust mass is calculated from the mass and metallicity of the cold 
gas in the galaxy, assuming that the dust-to-gas ratio is proportional 
to the metallicity. Since the gas mass and metallicity both evolve, 
so does the dust mass, and this evolution is fully taken into account 
in GRASIL. For simplicity, we assume that the size distribution of 
the dust grains and PAH molecules does not evolve, apart from the 
normalization. 

Once we have calculated the SEDs for the model galaxies, we 
compute luminosities in different observed bands (e.g. the optical 
B-band or the Spitzer 24/mi band) by convolving the SED with the 
filter+detector response function for that band. For computing the 
predicted fluxes from galaxies in a fixed observer-frame band, we 
redshift the SED before doing the convolution. 

The GRASIL code is quite CPU-intensive, requiring several 
minutes of CPU time per galaxy. Consequently, we are limited to 
running samples of a few thousand galaxies at each redshift. As a 
result, quantities such as luminosity functions and redshift distri- 
butions still show some small amount of noise, rather than being 
completely smooth curves, as can be seen in many of the figures in 
this paper. 

2.3 Choice of parameters in the GALFORM+GRAS I L model 

The combined GALFORM+GRASIL model has a significant num- 
ber of parameters, but this is inevitable given the very wide range 
of physical processes which are included. The parameters are con- 
strained by requiring the model predictions to reproduce a limited 
set of observational data - once this is done, there is rather little 
freedom in the choice of parameters. We have described above how 
the main paramet ers a re fixed, and more d etails can be found in 
ICole et al. I d2000h and lBaugh et al.\ {2005). For both of these pa- 
pers, large grids of GALFORM models were run with different pa- 
rameters, in order to decide which set of parameters gave the best 
overall fit to the set of calibrating observational data. These papers 
also show the effects of varying some of the main model parameters 
around their best-fit values. The parameters in the standard model 
for which we present results in this paper were chosen to reproduce 
the following properties for present-day galaxies: the luminosity 



functions in the B- and K-bands and at 60/im, the relations between 
gas mass and luminosity and metallicity and luminosity, the size- 
luminosity relation for galaxy disks, and the fraction of spheroidal 
galaxies. In addition, the model was required to reproduce the ob- 
served rest-frame far-UV (1500A) luminosity function at z = 3, 
and the obs erved sub-mm n umber counts and redshift distribution 
at 850/im teaughef a/.l2005h . The sub-mm number counts are the 
main factor driving the need to include a top-heavy IMF in bursts. 

The parameters fo r our standard model are exactly the same as 
in iBaughef a/.U2005h . which were chosen before Spitzer data be- 
came available. Since these parameters were not adjusted to match 
any data obtained with Spitzer, the predictions of our model in the 
Spitzer bands are genuine predictions. We could obviously have 
fine-tuned our parameters in order to match better the observational 
data we considered in this paper, but this would have conflicted 
with our main goal, which is to present predictions for a wide set of 
observable properties based on a single physical model in a series 
of papers. 

Since our assumption of a top-heavy IMF in bursts is a con- 
troversial one, we will also show some predictions from a variant 
model, which is identical to the standard model, except that we 
assume the same solar neighbourhood (Kennicutt) IMF in bursts 
and in disks. Comparing the predictions for the standard and vari- 
ant models then shows directly the effects of changing the IMF 
in bursts. We note that the variant model matches the present-day 
optical and near-IR luminosity functions almost as well as the stan- 
dard model, though it is a poorer fit to the local luminosity 
function for the brightest galaxies (see Fig. [5). The variant model 
underpredicts the 850/im counts by a factor of 10-30. 



3 NUMBER COUNTS 

We begin our comparison of the predictions of our galaxy forma- 
tion model against Spitzer data with the galaxy number counts. 
Fig. shows number counts in the four IRAC bands (3.6, 4.5, 5.8 
and 8.0 /im), and Fig. [2]does the same for the three MIPS bands 
(24, 70 and 160 fim). Each panel is split in two: the upper sub- 
panel plots the counts per logarithmic flux interval, dN/ d\nS v , 
while the lower sub-panel instead plots S v dN / d\n S v . The latter 
is designed to take out much of the trend with flux, in order to 
show more clearly the differences between the model and the on- 
servational data. In each case we plot three curves for our standard 
model: the solid blue line shows the total number counts includ- 
ing both extinction and emission by dust, the solid red line shows 
the contribution to this from galaxies currently forming stars in a 
burst, and the solid green line shows the contribution from all other 
galaxies (star-forming or not), which we denote as "quiescent". In 
Fig. [T] we also plot a dashed blue line which shows the predicted 
total counts if we ignore absorption and emission from interstellar 
dust (emission from dust in the envelopes of AGB stars is still in- 
cluded in the stellar contribution, however). In the MIPS bands, the 
predicted counts are negligible in the absence of interstellar dust, 
so we do not plot them in Fig. [2] In the lower sub-panels, we also 
show by a dashed magenta line the prediction from a variant model 
which assumes a normal (Kennicutt) IMF for all star formation, but 
is otherwise identical to our standard model (which has a top-heavy 
IMF in bursts). This variant model fits the local B- and K-band and 
60 /im luminosity functions about as well as our standard model, 
but dramatically underpredicts the 850 /im number counts. The ob- 
served number counts are shown by black symbols with error bars. 
Overall, the agreement between the predictions of our stan- 



8 Lacey et al. 




Figure 1. Galaxy differential number counts in the four IRAC bands. The curves show model predictions, while the symbols with error bars show observational 
data from Fazio et al. 1 2004) (with different symbols for data from different survey fields). Each panel is split in two: the upper sub-panel plots the counts 
as dN /d\nS v vs S v , while the lower sub-panel plots S u dN /d\nS u (in units mjy deg -2 ) on the same horizontal scale. The upper sub-panels show four 
different curves for our standard model - solid blue: total counts including dust extinction and emission; dashed blue: total counts excluding interstellar dust; 
solid red: ongoing bursts (including dust); solid green: quiescent galaxies (including dust). The lower sub-panels compare the total counts including dust for 
the standard model (solid blue line) with those for a variant model with a normal IMF for all stars (dashed magenta line). The vertical dashed line shows the 
estimated confusion limit for the model, (a) 3.6 /im. (b) 4.5 /im. (c) 5.8 /im. (d) 8.0 /Ltm. 



dard model and the observed counts is remarkably good, when one 
takes account of the fact that no parameters of the model were ad- 
justed to improve the fit to any data from Spitzer. Consider first the 
results for the IRAC bands, shown in Fig. [T] Here, the agreement 
of the model with observations seems best at 3.6 and 8.0 /im, and 
somewhat poorer at 5.8 /im. The model predicts somewhat too few 
objects at fainter fluxes in all of the IRAC bands. Comparing the red 
and green curves, we see that quiescent galaxies rather than bursts 
dominate the counts at all observed fluxes in all of the IRAC bands, 
but especially at the shorter wavelengths, consistent with the expec- 
tation that at 3.6 and 4.5 /im, we are seeing mostly light from old 
stellar populations. Comparing the solid and dashed blue lines, we 
see that the effects of dust are small at 3.6 and 4.5 /im, with a small 
amount of extinction at faint fluxes (and thus higher average red- 
shifts), but negligible extinction for brighter fluxes (and thus lower 
redshifts). On the other hand, dust has large effects at 8.0 /im, with 



dust emission (due to strong PAH features at A ~ 6 — 9/tm) be- 
coming very important at bright fluxes (which correspond to low 
average redshifts - see Fig. IAl! b) in the Appendix). The 8.0 /im 
counts thus are predicted to be dominated by dust emission from 
quiescently star-forming galaxies, except at the faintest fluxes. The 
counts at 5.8 /im show behaviour which is intermediate, with mild 
emission effects at bright fluxes and mild extinction at faint fluxes. 
Comparing the solid blue and dashed magenta lines, we see that the 
predicted number counts in the IRAC bands are almost the same 
whether or not we assume a top-heavy IMF in bursts, consistent 
with the counts being dominated by quiescent galaxies. 

Consider next the results for the MIPS bands, shown in Fig.[2] 
We again see remarkably good agreement of the standard model 
with the observational data. The agreement is especially good at 
faint fluxes (corresponding to higher redshifts). In particular, the 
model matches well the observed 24 /im counts at the "bump" 
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around fluxes S u ~ 0.1 — lmjy. Accurate modelling of the PAH 
emission features is obviously crucial for modelling the 24 fim 
number counts, since the PAH features dominate the flux in the 
24 pxa band as they are redshifted into the band at z > 0.5. 
On the other hand, the standard model overpredicts the number 
counts at bright fluxes (corresponding to low redshifts) in all three 
MIPS bands. The evolution at these wavelengths predicted by our 
ACDM-based model thus seems to be not quite as strong as indi- 
cated by observations. 

In the MIPS bands, emission from galaxies is completely 
dominated by dust, which is why no dashed blue lines are shown in 
Fig. [2] Comparing the red and green curves, we see that quiescent 
(but star-forming) galaxies tend to dominate the number counts in 
these bands at brighter fluxes, and bursts at fainter fluxes. This re- 
flects the increasing dominance of bursts in the mid- and far-IR 
luminosity function at higher redshifts. Comparing the solid blue 
and dashed magenta curves, we see that our standard model with a 
top-heavy IMF in bursts provides a significantly better overall fit to 
the observed 24 fim counts than the variant model with a normal 
IMF in bursts (although at the brightest fluxes, the variant model 
fits better). The faint number counts at 70 /jm also favour the top- 
heavy IMF model, while the number counts at 160 fim cover a 
smaller flux range, and do not usefully distinguish between the two 
variants of our model with different burst IMFs. 

We can use our model to predict the flux levels at which 
sources should become confused in the different Spitzer bands. 
We estimat e the confusion limit using the source density cri- 
terion (e.g. IVaisanen et al. I feOQll: iDole et al. 1 120031) : if the tele- 
scope has an FWHM beamwidth of Qfwhm, we define the ef- 
fective beam solid angle as ujbeam = (7r/(41n2)) Ofwhm = 
1.1W FWHM , and then define the confusion limited flux S CO nf 
to be such that N(> S conf ) = l/(N beam w beam ), where 
iV(> S) is the number per solid angle of sources brighter than 
flux S. We choose Nbeam = 20 for the number of beams 
per sour ce, which gives simi l ar results to more detailed analy- 
ses (e.g. IVaisanen et al. I [200 ll ; IDole et al"ll2004bl) . We use values 
of the beamsize 9fw hm = (1.66, 1.72, 1.88, 1.98) arcsec for 
the four IRAC bands dFazio et al. I l2004bl) and (5.6 16.7, 35.2) 
arcsec for the three MIPS bands {DojeefoT] [2003). Our stan- 
dard model then predicts confusion-limited fluxes of S con f = 
(0.62, 0.62, 0.69, 0.70)MJy in the (3.6, 4.5, 5.8, 8.0)^m IRAC 
bands, and Sconf = (0.072, 2.6, 43)mJy in the (24, 70, 160) Aim 
MIPS bands. Thes e confusion estimat es for the MIPS bands are 
similar to those of IDole et al. 1 l l2004bh . which were based on ex- 
trapolating from the observed counts. These values for the confu- 
sion limits are indicated in Figs.Q]and[2]by vertical dashed lines. 

Our galaxy evolution model does not compute the contribution 
of AGN to the IR luminosities of galaxies. On the other hand, the 
observed number counts to which we compare include both normal 
galaxies, in which the IR emission is powered by stellar popula- 
tions, and AGN, in which there is also IR emission from a dust 
torus, which is expected to be most prominent in the mid-IR. How- 
ever, multi-wavelength studies using optical, IR and X-ray data in- 
dicate that even at 24 pm, the fraction of so urces dominated at 
that w avelength by AGN is only 10-20% (e.g. iFranceschini et al. I 
l2005h . and the contribution of AGN-dominated sources in the other 
Spitzer bands is likely to be smaller. Therefore we should not make 
any serious error by comparing our model predictions directly with 
the total number counts, as we have done here. 
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Figure 2. Galaxy differential number counts in the three MIPS bands. The 
curves show model predictions while the symbols with error bars show ob- 
servational data. The meaning of the different mo del lines is the same a s 
in Fig. [TJ (a) 24 /im, with observational data fromlPapovich et al. I J2004I) . 
(b) 7 (im, with observati onal data fro mlDole et al. I J2004al) (filled sym- 
bols). lFraver et q/.l j2006al) (crosses), and lFraver et al. I j2006bl) (open sym- 
bols), ( c) 160 fim (bottom p anel), with observatio nal data from lDole et al. I 
1 2004a) (filled symbols) and Frayer et al. (2006a) (crosses). 
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4 EVOLUTION OF THE GALAXY LUMINOSITY 
FUNCTION 

While galaxy number counts provide interesting constraints on the- 
oretical models, it is more physically revealing to compare with 
galaxy luminosity functions, since these isolate behaviour at partic- 
ular redshifts, luminosities and rest-frame wavelengths. In the fol- 
lowing subsections, we compare our model predictions with recent 
estimates of luminosity function (LF) evolution based on Spitzer 
data. 



4.1 Evolution of the galaxy luminosity function at 3-8 pm 

We consider first the evolution of the luminosity function in the 
wavelength range covered by the IRAC bands, i.e. 3.6-8.0 (im. 
Fig. [3] shows what our standard model with a top-heavy IMF in 
bursts predicts for LF evolution at rest-frame wavelengths of 3.6 
and 8.0 /im for redshifts z = — 30- We see that at a rest-frame 
wavelength of 3.6 (im, the model LF hardly evolves at all over the 
whole redshift range z = — 3. This lack of evolution appears to 
be somewhat fortuitous. Galaxy luminosities at a rest-frame wave- 
length of 3.6 are dominated by the emission from moderately 
old stars, but the stellar mass function in the model evolves quite 
strongly over the range z — — 3 (as we show in Sj5). The weak 
evolution in the 3.6 pm LF results from a cancellation between a 
declining luminosity-to-stellar-mass ratio with increasing time and 
increasing stellar masses (see Figs. [T3l a) and (e)). On the other 
hand, at a rest-frame wavelength of 8.0 pm, the model LF becomes 
significantly brighter in going from z — to z = 3. Galaxy lu- 
minosities at a rest-frame wavelength of 8.0 pm are dominated by 
emission from dust heated by young stars, so this evolution reflects 
the increase in star formation activity with increasing redshift (see 
Fig. [Hb) in fg). 

In Fig. [4] we compare the model predictions for evolu- 
tion of the LF at 3.6 /im with observational estimates from 



3 In this figure, and in Figs.[4][5][8] and |10| the luminosities L„ are calcu- 
lated through the corresponding Spitzer passbands. 



iBabbedge etal] d2006h and iFranceschini et al. I |2006D The 
model predictions are given for redshifts z = 0, 0.5 and 1. For 
the observational data, the mean redshifts for the different redshift 
bins used do not exactly coincide with the model redshifts, so we 
plot them with the model output closest in redshiftQ. The observa- 
tional estimat es of the 3.6 LF re ly on the measured redshifts. 
In the case of IBabbedge et al\ |2006), these are mostly photomet- 
ric, usi ng optical and NIR ( including 3.6 and 4.5 (im) fluxes, while 
for the IFranceschini et al. I sample, about 50% of the redshifts are 
spectroscopic and the remainder photometric. In both samples, the 
measured 3.6 pm fluxes were k-corrected to estimate the rest-frame 
3.6 pm luminosities. 

We see from comparing the blue curve with the observational 
data in Fig.|4]that the 3.6 /an LF predicted by our standard model 
is in very good agreement with the observations. In particular, the 
observational data show very little evolution in the 3.6 pm LF over 
the redshift range z = — 1. Th e largest difference seen is at 
z = 1, where the Babbedg e et al. I data show a tail of objects to 
very high luminosities, which is not seen in the model p redictions. 
However, this tail is not seen in the IFranceschini "e t al. I data at the 
same redshift, and is also not present in the observational data at 
the lower re dshifts. More spectroscopic redshifts are needed for the 
Babbedge et al. sample to clarify whether this high-luminosity tail 
is real. Comparing the red, green and blue lines for the standard 
model shows that the model luminosity function is dominated by 
quiescent galaxies at low luminosity, but the contribution of bursts 
becomes comparable to that of quiescent galaxies at high luminosi- 
ties. We have not shown model LFs excluding dust extinction in 
this figure, since they are almost identical to the predictions includ- 
ing dust. The dashed magenta lines show the predicted LFs for the 

4 IBabbedge efoTlfeOOel) also compared their measured LFs at 3.6, 8.0 and 
24 pm with predictions from a preliminary version of the model described 
in this paper 

5 Specifically, f or z = 0, we compare with the z = 0.1 data from 
IBabbedge et al. L for z = 0.5 we compar e with the z = .5 data 
from [Babbedge et al.~\ and z = 0.3 data from Franceschini et al. , and for 
z = 1, we compare with t he z = 0.75 (op en symbols) and z = 1.25 
(filled symbols) da ta from IBabbedge et al. I and z = 1.15 data from 
IFranceschini et aT\ 
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Figure 4. Predicted evolution of the galaxy luminosity function at rest- 
frame 3.6 (im compared to observational data. The different panels show 
redshifts (a) z = 0, (b) z = 0.5 and (c) z = 1. The predictions for our stan- 
dard model are shown by the blue line, with the red and green lines showing 
the separate contributions from ongoing bursts and quiescent galaxies. The 
dashed magenta line shows the prediction for a variant model with a nor- 
mal IMF for all stars. The error bars on the model lines indicate the Poisson 
uncertainties due to the finite number of galaxies simulated. The black sym- 



Figure 5. Predicted evolution of the galaxy luminosity function at rest- 
frame 8.0 £tm compared to observational data. The different panels show 
redshifts (a) z = 0, (b) z = 1 and (c) z = 2. The coloured lines showing 
the model predictions have the same meaning as in Fig. [4] T he black sym- 



bols with error bars show observational data from Babbedge et al. 
(open circles for z = and 0. 7, triangles for z = 1.2), lHuang et al. 
(filled circles for z = 0) and lCaputi et al. I fe007l) (filled circles for 
and 2). The observed LFs are for normal galaxies and exclude AGN 
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variant model with a normal IMF in bursts. We see that these differ 
only slightly from our standard model, but are a somewhat poorer 
fit to the observational data at higher luminosities. 

In Fig.[5]we show a similar comparison for the LF evolution at 
a rest-frame wavelength of 8 (jm. The model predictions are given 
for redshifts z = 0, 1 and 2 , and are compared w ith observational 
estima tes by iHuang et al] < |2007|) (for z ~ 0). iBabbedge et al. I 
b006h (for z ~ and z ~ 1) and lCaputi etal] d2007t) (for z ~ 1 
and z ~ 2). These papers all classified objects in their samples 
as either galaxies or AGN, and then computed separate LFs for 
the two types of objects 0. Our model does not make any predic- 
tions for AGN, so we compare our model predictions with the ob- 
served LFs for objects classified as galaxies only . We see that for 
red shifts around z = 1, the observed LFs from Babbed ge et al. I 
an dlCaputi et al. I are in very poor agreement with each other, with 
the lCaputi et al. I LF being around 10 times higher in number den- 
sity at the same luminosity. This differerence presumably results 
from some combination of: (a) different met hods of classifying ob- 
jects as galaxies or AGN (Babb edge et al. I used only optical and 
IR fluxes to do this, while ICaputi et al. also used X-ray data); (b) 
different photometric redshift estimators; and (c) different meth- 
ods for k-correcting luminosities to a rest-fram e wavelength of 8 
^m. There are s maller differences between the IHuang et al I and 
Babbedge et al. LFs at z ~ 0. Futher observational investigation 
appears to be necessary to resolve th ese issues. Our sta ndard model 
is in reasonable agreem ent with the [Babbedge et al. I observed LF 
at z ~ 0, and with th e ICaputi et al. I ob served LFs at z ~ 1 and 
2 ~ 2, but not wi th the Babbedge et al. I observed LF at z ~ 1. The 
comparison with ICaputi et al. favours our standard model with a 
top-heavy IMF in starbursts over the variant model with a normal 
IMF. 



4.2 Evolution of the galaxy luminosity function at 12-24 /im 

In this subsection, we consider the evolution of the galaxy lumi- 
nosity function at mid-IR wavelengths, and compare with data ob- 
tained using mainly the MIPS 24 /im band. 

Fig. [6] shows what our standard model with a top-heavy IMF 
in bursts predicts for the evolution of the galaxy LF at rest-frame 
wavelengths of 15 and 24 fim for redshifts z — — 3 01 At rest- 
frame wavelengths of 15 and 24 /im, galaxy luminosities are typi- 
cally dominated by the continuum emission from warm dust grains 
heated by young stars (although PAH emission is also significant 
at some nearby wavelengths). Fig. [6] shows strong evolution in the 
model LFs over the redshift range z = — 3 at both wavelengths, 
reflecting both the increase in star formation activity with increas- 
ing redshift (see Fig. [l~3l b)) and the increasing dominance of the 
burst mode of star formation, for which the top-heavy IMF fur- 
ther boosts the mid- and far-IR luminosities compared to a normal 
IMF. Comparing Fig. [6] with Fig.Oa), we also see a difference in 
the shape of the bright end of the LF: at 3.6 fim, where the LF is 
dominated by emission from stars, the bright end cuts off roughly 



6 Note that a variety of criteria have been used for classifying observed IR 
sources as AGN or normal galaxies, and these do not all give equivalent 
results. Even if an object is classified as an AGN, it is also not clear that in 
all cases the AGN luminosity dominates over that of the host galaxy in all 
Spitzer bands 

7 In this figure, and in Figs.[7]and[8] the 24 /im luminosities are calculated 
through the corresponding MIPS passband, while the 15 fim luminosities 
are calculated through a top-hat filter with a fractional width of 10% in 
wavelength. 
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Figure 8. Predicted evolution of the galaxy luminosity fun ction at rest- 
frame wavelength 24 pm compared to obse rvational data fro m Shupe et al. 
1 19981) (at 2 = 0, open symbols) and from lBabbedge et al. 1 120061) (for the 
same redshifts as in Fig. [4). The meaning of the curves showing the model 
predictions is the same as in Fig. [4] (a) z = 0, (b) z = 0.5 and (c) 2 = 1. 
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Figure 6. Predicted evolution of the galaxy luminosity function in our standard model at rest-frame wavelengths(a) 15 pm (left) and (b) 24 pm (right) for 
redshifts z = 0, 0.5, 1, 1.5, 2 and 3, as shown in the key. 



exponentially, while at 15 and 24 /im, where the LF is dominated 
by emission from warm dust, the bright end declines more gradu- 
ally, roughly as a power-law. This difference reflects the difference 
in shape of the galaxy stellar mass function (GSMF) and the galaxy 
star formation rate distribution (GSFRD). The GSMF shows an 
exponential-like cutoff at high masses, while the GSFRD shows a 
more gradual cutoff at high SFRs because of starbursts triggered by 
galaxy mergers (see Figs.|13fa) and (b) in Sj5). This difference was 
noticed earlier by observers comparing opt ical and far-IR LFs of 
galaxies, but its or igin was not understood dLawrence et al. Ifl986l : 
ISoifere? aZ. Ill 987bl) . 

In Fig. [7] we compare the model LFs at rest-frame wave- 
lengths 12 and 15 fim with observational estimates. For z = 0, 
we pl ot th e observational es timates from Soifer & Neugcbaucr 
dl99ll) and lRush etaT\ dl993h . based on IRAS 12 /im data (with 
AGN removed ). For z = 0.5 — 1 and z = 1.5 — 2.5, we 
plot the data o Floc'h et al. I d2005l) and fPerez -Gonzalez et al. I 
(2005) respectively, which were o btained from gal axy samples 
selected on Spitzer 24 /im flux. iLe Floc'h et al. I k-corrected 
their measured 24 pm fluxe s to 15 pm rest-frame luminosities, 
while | Perez-Gonzalez et aT] k-corrected to 12 pm rest-framfl 
llx Floc'h et al. I obtained most of th eir redshifts from photo metric 
redshifts based on optical data, while lPerez-Gonzalez et aZ.l used a 
new photometric redshift technique based on fitting empirical SEDs 
to all of the available broad-band data from the far-UV to 24 pm, 
and also removed "extreme" AGN from their observed LF. Note 
that the redshifts for the observed LFs do not exactly coincide with 
model redshifts in all cases, but are close. 

We see from comparing the blue line to the observational dat- 
apoints in Fig.|7]that our standard model with a top-heavy IMF in 
bursts fits the observations remarkably well up to z — 2. In partic- 
ular, the model matches the strong evolution in the mid-IR LF seen 



8 The exact passband used for the model LF in each panel depends on 
which observational data we are comparing with. For z = 0, we use the 
IRAS 12 pm passband; at z = 0.5 and z = 1 we use a top-hat passband 
centred at 15 pm; and at z = 1.5, 2 and 2.5, we use a top-hat passband 
centred at 12 pm (both top-hat passbands having fractional width 10% in 
wavelength). 



( 1998) (for 
20061) (for 



in the observational data. The model falls below the observational 
data at z = 2.5, but here both the photometric redshifts and the 
k-corrections are probably the most uncertain. The standard model 
also does not provide a perfect fit to the z = data, predicting 
somewhat too many very bright galaxies and somewhat too few 
very faint galaxies (though the latter discrepancy might be affected 
by local galaxy clustering in the IRAS data). Comparing the red, 
green and blue lines for the standard model in the figure, we see 
that the bright end of the 12 or 15 /im LF is dominated by bursts 
at all redshifts. The figure also shows by a dashed magenta line the 
predictions for the variant model with a normal IMF in bursts. This 
latter model predicts much less evolution in the bright end of the LF 
than is observed. This comparison thus strongly favours the model 
with the top-heavy IMF in bursts. 

Finally, in Fig.(8] we carry out a similar comparison of the evo- 
lution of predicted and observed LFs at a rest-frame wavelength 
of 24 pm over the redshift range z = — 1, in this case com- 
paring with observational estimates fr omlShupe et al. 
z — 0), based on IRAS data, and from Babb edge et al. _ 
z = — l), bas ed on Spitzer datfl The galaxy redshifts for the 
Babbedge et al. data were obtained in the same way as for the 3.6 
pm LFs shown in Fig. [4] and the luminosities were k-corrected 
from observer-frame 2 4 pm to rest-frame 24 pm. The LF plotted 
from [Babbedge et al. I is that for normal galaxies, with AGN ex- 
cluded. 

The conclusions from comparing the model with the 24 
LFs are similar to those from the comparison with the 12 and 15 
/im LFs. The data favour our standard model over the variant with 
a normal IMF in bursts (except possibly for z — 0.5), as the latter 
predicts too little evolution at the bright end. At z — 0, the model 
fits the 24 pm data rather better than for the corresponding com- 
parison at 12 pm. On the other hand, at z — 0.5 and z = 1, the 
model LF is a somewhat worse fit to the observational data at 24 
pm than at 15 pm. These differences between the 12/15 and 24 
pm comparisons might result from the different photometric red- 
shifts and k-corrections used in the observational samples in the 



■' The model luminosities are all computed through the Spitzer 24 pm pass- 
band. 
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Figure 7. Predicted evolution of the galaxy luminosity function at rest-frame wavelength 12 or 15 /an compared to observational data. The different panels 
show redshifts: (a) z = 0, (b) z = 0.5, (c) z = 1 , (d) z = 1.5, (e) z = 2 and (f) 2 = 2.5. The meaning of the curves showing t he model predictions is the 
same as in F ig. [4] In panel (a), the predictions at 12/xm are compared to observational determinations from lSoifer & Neugebauer U1991I) (op en symbols) and 
iRush et al. (1993) (filled symbols) based on IRAS data. In panels (b) and (c), the predictions at 15/^ m are compared to observatio nal data from lLe Floc'h et al. I 
boost) . In panels (d), (e) and (f), the predictions at 12/im are compared to observational data from lPerez-Gonzalez et al. (2005). 
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Figure 9. The predicted galaxy luminosity function at 60 /an compared 
to observational data from IRAS. The meaning of the different lines is 
the same as in Fig. [4| The blac k symbols show observationa l data from 
ISaunders et al. I )l990l) (crosses) , ISoi fer & Neugebauer] jl99lh (open cir- 
cles). and lTakeuchi e? 0/71^20031) (filled circles). 



two cases. Alternatively, they might result from problems in mod- 
elling the dust SEDs in the complex mid-IR range. 



4.3 Evolution of the galaxy luminosity function at 70-160 /im 

We now briefly consider the evolution of the luminosity function 
in the far-IR. The far-IR is the wavelength range where most of 
the luminosity from dust in normal galaxies is emitted. The lo- 
cal 60 )im luminosity function was very well measured by sur- 
veys with IRAS, and so is commonly used as a starting point or 
benchmark for modelling the evolution of the galaxy population in 
the far-IR. We therefore present in Fig. [9] the model prediction for 
the 60 fj.m luminosit y function at z = , compared with obser- 
vation al da ta from ISaunders et al. l dl990h . lSoifer & Neugebauer 
199 lh and iTakeuchi et al. I J2003L As discussed in Bau gh et al. I 



2005), the local 60 /an luminosity function was used as one of the 



primary constraints in fixing the parameters of our galaxy forma- 
tion model, and the figure shows that our standard model provides 
a good match to the data. The variant model with a normal IMF in 
bursts underpredicts the abundance of the brightest 60 fim galaxies. 

In Fig. [10] we show the model predictions for the evolution of 
the luminosity function in the two longer wavelength MIPS bands, 
at rest-frame wavelengths of 70 and 160 /im, from z — to z = 3. 
At 70 /Urn, the luminosity function at high luminosities is predicted 
to brighten by about a factor 10 going from z = to z = 2. This 
is about a factor 2 less than the brightening predicted in the mid-IR 
at 15 p,m (compare to Fig.[6jl, but nearly a factor 2 more evolution 
than is predicted at 160 /im. These differences between the amount 
of evolution seen at different IR wavelengths reflect evolution in 
the shapes of the SEDs of the galaxies responsible for the bulk of 
the IR emission. No observational estimates of the evolution of the 
luminosity function at 70 and 160 p,m have yet been published, 
but they are expected to be forthcoming from ongoing surveys with 
Spitzer. 



Figure 11. Predicted evolution of the total mid+far-IR (8-1000 fim) galaxy 
luminosity function for our standard model, for redshifts z = 0, 1, 2, 3, 4 
and 6, as shown in the key. 



4.4 Evolution of the total mid+far-IR luminosity function 

The total mid+far IR luminosity of a galaxy, Lir, integrated over 
the whole wavelength range 8-1000 fim, is a very good approxi- 
mation to the total luminosity emitted by interstellar dust grains in 
all galaxies except those with very small dust contents. In galax- 
ies with significant star formation, Lir is mostly powered by dust 
heated by young stars, and so provides a quantitative indicator of 
the amount of dust-obscured star formation which is independent 
of the shape of the IR SED (though still subject to uncertainties 
about the IMF). The evolution of the luminosity function in Lir 
is therefore a very interesting quantity to compare between models 
and observations. We show in Fig.[TT]what our standard model pre- 
dicts for the evolution of the IR LF over the range 2 = — 6. We 
see that the model predicts substantial evolution in this LF, with the 
high luminosity end brightening by a factor ~ 10 from z — to 
z = 2, followed by a "plateau" from z = 2 to z = 4, and a decline 
from z = 4 to z — 6. 

In Fig. [J_2] we compare our model predictions with existing 
observational estimates of the total IR LF for 2 = — 2. These 
observational estimates are only robust for z = 0, where they are 
based on IRAS measurements covering the wavelength range 12- 
100 /im. At all of the higher redshifts plotted, the observational 
estimates are based on measurements of the mid-IR luminosity de- 
rived from Spitzer 24 /im fluxes, converted to total IR luminosi- 
ties by assuming SED shapes for the mid- to far-IR emission. The 
bolometric correction from the observed mid-IR luminosity to the 
inferred total IR luminosity is typically a factor ~ 10, and is sig- 
nificantly uncertain. Therefore, the most robust way to compare 
the models with the observations is to compare them at the mid- 
IR wavelengths where the measurements are actually made, as we 
have done in $4. 1 1 and £|4.2| Nonetheless, if we take the observa- 
tional determinations at face value, then we see that observed evo- 
lution of the total IR LF agrees remarkably well with the predic- 
tions of our standard model with a top-heavy IMF. On the other 
hand, the variant model with a normal IMF predicts far too few 



16 Lacey et al. 




log^Vh-* Lb) log^/h"* L G ) 

Figure 10. Predicted evolution of the galaxy luminosity function in our standard model (including dust) at rest-frame wavelengths (a) 70 (im and (b) 160 /im, 
for redshifts 2 = 0,0.5,1,1.5,2 and 3, as shown in the key. 



high Lie. galaxies at higher z, and is strongly disfavoured by the 
existing data. 



5 INFERRING STELLAR MASSES AND STAR 
FORMATION RATES FROM Spitzer DATA 

In this section, we consider what the models imply about how well 
we can infer the stellar masses and star formation rates (SFRs) in 
galaxies from measurements of rest-frame IR luminosities. The top 
two panels of Fig.|13|show the predicted galaxy stellar mass func- 
tion (GSMF, left panel) and galaxy star formation rate distribution 
(GSFRD, right panel), for redshifts z = — 6. We see that the pre- 
dicted stellar mass function shows dramatic evolution over this red- 
shift range, with a monotonic decline in the number of high-mass 
galaxies with increasing redshift. On the other hand, the SFR distri- 
bution shows much less dramatic evolution over this redshift range, 
with a mild increase in the number of high-SFR objects up to z ~ 3, 
followed by a decline above that. The lower four panels in Fig. [T3l 
show the relation in the models between stellar masses and SFRs 
and rest-frame luminosities at different IR wavelengths. (Note that 
in all cases, luminosities are measured in units of the bolometric 
solar luminosity.) The middle and bottom left panels respectively 
show the mean ratio of luminosity in the rest-frame K (2.2/im) or 
3.6 fj,m bands to stellar mass as a function of stellar mass. The 
middle and bottom right panels respectively show the mean ratio 
of total mid+far-IR (8 — lOOO^im) or rest-frame 15 /im luminosity 
to SFR as a function of SFR. (The mean L/M, or L/SFR ratios 
plotted are computed by dividing the total luminosity by the total 
mass or SFR, in each bin of mass or SFR.) 

The near-IR luminosity is often used as a tracer of stellar mass. 
The left panels of Fig,[T3lshow that the L/M, ratio varies strongly 
with redshift, reflecting the difference in the ages of the stellar pop- 
ulations. At higher redshifts it also shows a significant dependence 
on stellar mass, presumably reflecting a trend of age with mass. 
However, the variation of mean L/M, with redshift is seen to be 
much smaller at 2.2 /im than at 3.6 /im, implying that the rest- 
frame K-band light should provide a more robust estimator of stel- 
lar mass than the light at longer wavelengths. The differences be- 
tween L/M, values at 2.2 /im and 3.6 reflect the larger contri- 



bution from AGB compared to RGB stars at the longer wavelength. 
AGB stars have higher masses and younger ages than RGB stars, 
and so are more sensitive to star formation at recent epochs. The 
scatter in L/M, at a given mass is also found in the models to in- 
crease with redshift. In the K-band, it increases from ~ 40% at 
z ~ to a factor ~ 3 at z ~ 6. The large scatter at high redshifts 
results in part from having two different IMFs. 

The luminosity in the mid- and far-IR is widely used as a 
tracer of dust-obscured star formation (although in galaxies with 
very low star formation rates, the dust heating can be dominated 
by older stars). The total mid+far-IR (rest-frame 8-1000 fj,m) lumi- 
nosity is expected to provide a more robust tracer of star formation 
than the luminosity at any single IR wavelength, since the shape of 
the SED of dust emission depends on the dust temperature distri- 
bution (as well as on the dust grain properties). This is borne out 
by our model predictions. The middle right panel of Fig.|13|shows 
that the Lir/SFR ratio depends weakly on both SFR and red- 
shift. This behaviour results mostly from having different IMFs in 
the model in quiescent and bursting galaxies, with the fractional 
contribution of the bursts increasing both with SFR and with red- 
shift. If we look at quiescent and bursting galaxies separately, we 
find roughly constant ratios Lir/SFR fa 6 x 10 9 h~ 1 L Q /M e 
and Lir/SFR = 2 x 1O 1o /i _1 L /M Q respectively, for galaxies 
where Lir is powered mostly by young stars. However, there is 
also a trend at lower redshift for Lir/SFR to be larger at lower 
SFR - this reflects the larger fraction of dust heating from older 
stars in galaxies with lower SFRs, which more than compensates 
for the lower average dust obscuration in these galaxies. The lower 
right panel of Fig.|13|shows that the L/SFR ratio in the mid-IR 
(in this case at 15 (j,m in the rest-frame) shows more variation with 
SFR and redshift than the ratio for the total IR luminosity. This re- 
flects the variation in the mid- to far-IR SED shapes in the model. 
The scatter in the L/SFR ratio is roughly a factor 2 around the 
average relation for the total IR luminosity, but is larger for the 15 
fim luminosity. 

The results of this section illustrate why it is not straightfor- 
ward to compare theoretical predictions for the evolution of the 
galaxy stellar mass function and star formation rate distribution (or 
even the stellar mass and star formation rate densities) with obser- 
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Figure 12. Predicted evolution of the total mid+far IR (8-1000/^m) galaxy LF compared to observational data. T he different panels show redshi fts (a) z = , 
(b) z = 0.5, (c) z = 1 and (d) 2 = 2. For z = 0, we compare with observational data from Sanders et at. 1 2003) (filled symbols) and Takeuchi et al. ( 2003) 
(open symbols, convert ing his 60 (im LF to a total IR LF assuming a constant c onversion factor, Lif j/vL v (60iim) = 2.5). We compare with data from 
iLe Floc'h e t al. (2005) for z = 0.5 and 2 = 1 (filled and open symbols), and with lCaputi et al. 1 120 07) for 2 = 1 and 2 = 2 (crosses). 



vational estimates. In addition to assumptions about galaxy star for- 
mation histories and metallicities (for stellar mass estimates), and 
about the SED shapes for dust emission (for SFR estimates from 
IR and sub-mm data), observational estimates all rest on some as- 
sumed form for the IMF. If the IMF assumed in the observational 
analysis is different from the true IMF, the observational estimates 
for stellar masses and SFRs can be wrong by large factors. If the 
IMFs differ only below IMq, then one can apply a simple rescal- 
ing to relate stellar mass and SFR estimates for different IMFs. 
However, if our current galaxy formation model is correct, stars 
form with different IMFs in quiescent disks and in merger-driven 
bursts, and so no observational estimate based on assuming a sin- 
gle IMF can give the correct GSMFs and GSFRDs, nor the correct 
stellar mass and SFR densities. A direct comparison of the GSMF 
and GSFRD evolution predicted by our model with observational 
estimates is therefore not meaningful. Instead, the comparison be- 
tween models and observations must be made via directly observ- 
able (rather than inferred) quantities, such as the K-band luminosi- 



ties to constrain stellar masses, and the total IR luminosities to con- 
strain SFRs. 



6 CONCLUSIONS 

We have computed predictions for the evolution of the galaxy pop- 
ulation at infrared wavelengths using a detailed model of hierar- 
chical galaxy formation and of the reprocessing of starlight by 
dust, and compared these predictions with observational data from 
the Spitzer Space Telescope. We calculated galaxy formation in 
the framework of the ACDM model using the GALFORM semi- 
analytical model, which includes physical treatments of the hier- 
archical assembly of dark matter halos, shock-heating and cool- 
ing of gas, star formation, feedback from supernova explosions 
and photo-ionization of the IGM, galaxy mergers and chemical en- 
richment. We computed the IR luminosities and SEDs of galaxies 
using the GRASIL multi-wavelength spectrophotometric model, 
which computes the luminosities of the stellar populations in galax- 
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Figure 13. Model predictions for properties related to stellar masses (left column) and star formation rates (right column), for redshifts z = 0, 1, 2, 3, 4, and 
6: (a) galaxy stellar mass function (GSMF); (b) galaxy star formation rate distribution (GSFRD); (c) mean ratio of rest-frame K-band luminosity to stellar 
mass, as a function of stellar mass; (d) mean ratio of total mid+far IR luminosity to SFR, as a function of SFR; (e) mean ratio of rest-frame 3.6 fim luminosity 
to stellar mass, as a function of stellar mass; (f) mean ratio of rest-frame 15 fim luminosity to SFR, as a function of SFR. (The 15 (im luminosity is here 
calculated through top-hat filter with a fractional wavelength width of 10%.) 
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ies, and then the reprocessing of this radiation by dust, including 
radiative transfer through a two-phase dust medium, and a self- 
consistent calculation of the distribution of grain temperatures in 
each galaxy based on a local balance between heating and cool- 
ing. The GRASIL model includes a treatment of the emission from 
PAH molecules, which is essential for understanding the mid-IR 
emission from galaxies. 

Our galaxy formation model incorporates two different IMFs: 
quiescent star formation in galaxy disks occurs with a normal so- 
lar neighbourhood IMF, but star formation in bursts triggered by 
galaxy mer gers happens with a top-heavy x = IMF. In a previ- 
ous paper (Baug hef al. 120051) . we found that the top-heavy IMF in 
bursts was required in order that the model reproduces the observed 
number counts of the faint sub-mm galaxies detected at 850 /im, 
which are typically ultra-luminous starbursts at z ~ 2, with total IR 
luminosities Lir ~ 10 12 — 1O 13 L0. This conclusion was arrived 
at following a search of a large grid of model parameters, with the 
imposition of a v ariety of deta i led ob servational constraints. The 
parameters in the iBaugh et al. I j2005h model were chosen before 
the publication of any results from Spitzer, without reference to any 
IR data apart from the local 60 /jm luminosity function and the 850 
/jm galaxy counts. We have kept the same parameter values in the 
present paper, in order to test what the same model predicts at other 
wavelengths and other redshifts. By doing this, we hope to address 
the criticism made of many semi-analytical models that they have 
no predictive power, because their parameters are always adjusted 
to match the observational data being analysed at that instant. 

We first compared the predictions from our model with the 
galaxy number counts measured in all 7 Spitzer bands, from 3.6 to 
160 fim. We found broad agreement between the model and the 
observations. In the 4 IRAC bands (3.6-8.0 /im), where the counts 
are mostly dominated by emission from older stellar populations, 
we found that the predicted counts were insensitive to whether we 
had a top-heavy or normal IMF in bursts. On the other hand, in 
the MIPS bands (24-160 /im), where the counts are dominated by 
emission from dust in star-forming galaxies, the predicted counts 
are more sensitive to the choice of IMF, and the counts are fit better 
by the model with a top-heavy IMF. We next investigated the evo- 
lution of the galaxy luminosity function at IR wavelengths, where 
several groups have now used Spitzer data to try to measure the 
evolution of the galaxy luminosity function over the redshift range 
z ~ — 2, at rest-frame wavelengths from 3.6 to 24 fim. 

Our model predicts that at mid- and far-IR rest-frame wave- 
lengths, the luminosity function evolution is very sensitive to the 
choice of IMF in bursts. We found that our standard model with a 
top-heavy IMF in bursts fits the measured evolution of the mid-IR 
luminosity function remarkably well (when allowance is made for 
complexity of predicting dust emission in the mid-IR), without any 
adjustment of the parameters. On the other hand, a model with a 
normal IMF in bursts predicts far too little evolution in the mid-IR 
luminosity function compared to what is observed. We made a sim- 
ilar comparison with the evolution of the total IR luminosity func- 
tion, where in the case of the observations, the total IR luminosities 
at high redshifts have been inferred from the 24 /im fluxes by fit- 
ting SEDs, and reached the same conclusion. The evolution of the 
galaxy luminosity function in the mid-IR found by Spitzer thus sup- 
ports our original conclusion about the need for a top-heavy IMF 
in bursts, which was based only on the sub-mm counts. This con- 
clusion will be further tested by ongoing Spitzer surveys at longer 
wavelengths. To assist this, we have also presented predictions for 
the evolution of the luminosity function in the Spitzer 70/im and 
160/im bands. 



We have also presented predictions for the evolution of the 
stellar mass function and star formation rate distribution of galax- 
ies. We investigated how the L/M* and L/SFR ratios varied with 
galaxy mass, SFR and redshift in different IR wavelength ranges, 
and considered the implications for observational estimates of stel- 
lar masses and SFRs from IR observations. Even in the near-IR, the 
predicted variations in L/ M* with mass and redshift can be surpris- 
ingly large. The variations mL/M, are much larger at a rest-frame 
wavelength of 3.6 (im than at 2.2 /im, implying that the 2.2 fim 
luminosity is a more robust tracer of stellar mass. 

Finally, we have presented in an Appendix the predictions of 
our model for the redshift distributions of galaxies selected at dif- 
ferent IR fluxes in the Spitzer bands. 

One significant limitation of our model is that it does not in- 
clude the effects of AGN. Two effects are relevant here. The first is 
feedback from AGN on galaxy formation. In several recent galaxy 
formation models, AGN feedback is invoked to prevent the forma- 
tion of too many massive galaxies at the present day. In the model 
presented here, we instead posit feedback from supernova-driven 
galactic superwinds, which perform a similar role to AGN feed- 
back in suppressing the formation of very massive galaxies. Both 
the superwind and AGN feedback models include free parameters 
which are tuned to give a match to the present-day optical galaxy 
luminosity function. However, the redshift dependence of the feed- 
back will be different between our superwind model and the various 
AGN feedback models, so in general they will all predict different 
evolution of the galaxy population with redshift. We will investi- 
gate galaxy evolution in the IR in a model with AGN feedback in 
a future paper. The second effect of AGN which we have not in- 
cluded is the emission from AGN and their associated dust tori. In 
order to compensate for this, we have wherever possible compared 
our model predictions with observations from which the AGN con- 
tribution has been subtracted out. This was possible for most of our 
comparisons of luminosity function evolution. This was not possi- 
ble for the number counts comparisons, but in this case the contri- 
bution from AGN is thought (based on observations) to be a small 
fraction of the total over the flux range explored by Spitzer, even in 
the mid-IR where the dust tori are the most prominent. We therefore 
believe that emission from AGN does not seriously affect our con- 
clusions about the IR evolution of star-forming galaxies. We hope 
to include AGN emission directly into our models in the future. 

We have thus shown that Spitzer data provide a stringent test of 
galaxy formation theory, by probing galaxy evolution, constraining 
star formation rates and the role of dust to z ~ 2. We find that 
an ab initio ACDM model gives an acceptable fit to the Spitzer 
data provided that ~ 10% of the stars in galaxies today formed 
in bursts of star formation with a top-heavy IMF. Future facilities 
like Herschel, SPICA, JWST and ALMA will continue to exploit the 
valuable information on galaxy formation contained in the IR part 
of the electromagnetic spectrum. 
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APPENDIX A: REDSHIFT DISTRIBUTIONS 

In this Appendix, we present some predictions from our standard 
model for the redshift distributions of galaxies selected at different 
fluxes in the Spitzer bands. This is principally for completeness, 
to assist in interpreting data from current surveys, and to assist in 
planning future surveys based on Spitzer data. The set of redshift 
distributions at all observed fluxes in principle contains equivalent 
information to that in the luminosity functions at different wave- 
lengths and redshifts. However, comparing models with observa- 
tions via luminosity functions is more physically transparent than 
making the comparison via redshift distributions, which is why we 
have presented our results on luminosity functions in the main part 
of the paper, and why we make only a limited direct comparison 
with observed redshift distributions in this Appendix. In addition, 
if one only compares the predicted and observed redshift distribu- 
tions for galaxies above a single flux limit (e.g. the flux limit of 
a survey), this has less information than comparing the luminosity 
functions at different redshifts. 

We first show in Fig. lAll how the median redshift, and the 10- 
90 percentile range, are predicted to change with flux for galax- 
ies selected in one of the four Spitzer bands 3.6, 8.0, 24 or 70 
|im. While at most wavelengths the median redshift is predicted 
to increase smoothly and monotonically with decreasing flux, this 
is not true at 24 pm, where there is a bump around S v ~ 100/iJy. 
The structure seen for the 24 pm band as compared to the other 
wavelengths results from different PAH emission features moving 
through the band with increasing redshift. 

In Fig. lA2l we show the predictions from our standard model 
for the redshift distributions of galaxies in the four IRAC bands. 
For each band, we show the redshift distribution for galaxies se- 
lected to be brighter than S„ > 10/iJy in that band. The flux limit 
S v > 10/iJy has been chose n to match that in the obse rved deep 
sample selected at 3.6/im by iFranceschini et al. 1 d2006h . In each 
panel, the blue curve shows the predicted dN/dz for all galax- 
ies, normalized to unit area under the curve, and the red and green 
curves show the separate contributions of bursting and quiescent 
galaxies to the total. For 3.6/i m, the black line shows the ob- 
served redshift distribution from lFranceschini et al. 1 12006), which 
has also been normalized to unit area under the curve. We see that 
the observed redshift distribution peaks at a slightly higher redshift 
than in the model. However, the luminosity function evolution de- 
rived from this same sample is in reas onable agreement with the 
model, as was already shown in Fig. l4l IFranceschini et al. I ( 2006) 
note that the peak seen in their data at z ~ 0.8 is partly contributed 
by large-scale structures in the CDFS field. 

In Fig. lA3l we show predicted redshift distributions for galax- 
ies selected to be at a set of different fluxes in the four IRAC bands. 
The curves for the different fluxes are all normalized to have unit 
area as before, but in this figure the galaxies are selected to be at a 
particular flux, rather than being brighter than a certain flux. As one 
would expect, the typical redshift increases as the flux decreases. 

Figs. IA4I and IA5I show for the three MIPS bands the equiv- 
alent of Figs. |A2] and [A3] for the IRAC bands. In Fig. |A4] we 
show the predicted redshift distributions for galaxies brighter than 
a particular flux, where this flux limit is taken to be 83 /iJy at 24 
pm, 10 mjy at 70 pm and 100 mjy at 160 pm. The flux limit 
at 24 /im has been chosen to match that u sed in the deep obser- 
vation al sa mples oflLe Floc'h et al. I (2005), Perez-Gonzalez et al. I 
(2005) and ICaputi et al. I (12006a). while the flux limits at 70 
and 160 /im have been chosen to be roughly 3 times brighter than 
the source confusion limits in these bands. We see in Fig. I A5 1 that 
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Figure Al. Model predictions for the median redshift as a function of flux in four Spitzer bands, (a) 3.6 fim, 8.0 fim, (c) 24 fim, (d) 70 /im. In each panel, 
the median redshift for galaxies at each flux is shown by a solid line, and the 10- and 90-percentiles are shown by dashed lines. 



the redshift distributions at 24 fim show much more structure than 
at other wavelengths. This results from different PAH emission fea- 
tures moving through the 24 fim band with changing redshift. 

In Fig. IA4f a). we compare the predicted redshift dis- 
tribution at 24 fim with observational determin ations from 
IPerez-Gonzalez et al. I d2005l) (dashed black line) and Cap uti et al. I 
(2006a) (solid black line). The observed distributions have been 
separately normalized to unit area under the curve, as for the 
model distribution. Both observed distributions are based pri- 
marily on photometri c redshifts, but the photometric redshifts of 
Caput i et al. (2006a) are lik ely to be more accurate than those of 
Perez-Gonzalez et al. I ( 120051) . since the for mer are based on deeper 
optical and K-band data than the latter. (Perez-Go nzalez et al. I 
found optical counterparts with Bab < 24.7 or Rab < 23.7 for 
~ 70% of their S v (24pm) > 83^Jy sources, but relied on IRAC 
fluxes in deriving photo- z's for the rem aining ~ 30% of their sam- 
ple. On the other hand, ICaputi et al. I found K-band counterparts 
with K(Vega) < 21.5 for 95% of their S„(24^m) > 80/iJy 
sample, and derived photo-z's for essentially all of these sources 
using optical and K-band data a lone). Both observed distributions 
are similar, but the lCaputi et al. I distribution shows more structure. 



This is a combination of the effects of more accurate photometric 
redshifts but also a 9 times smaller survey are a, which mean s that 
fluctuations due to galaxy clustering are larger. ICaputi et al\ argue 
that the separate peaks at z ~ 0.7 and 1.1 result from large-scale 
structure, but that the bump at z ~ 1.9 results from PAH emis- 
sion features entering the observed 24 band. We see that the 
model also predicts peaks in the redshift distribution at z ~ 0.3, 
z ~ 1 and z ~ 2, which can be explained by different PAH fea- 
tures moving through the 24 fim band, although the z ~ 2 peak is 
more prominent than is seen in the observational data. Overall, the 
model redshift distribution at this flux limit is too skewed to high 
redshift compared to the observations, predicting too few galaxies 
at z ~ 0.5 — 1, and too many in the peak at z ~ 2. 

We investigate further this apparent discrepancy in the 24 fim 
redshift distribution in Fig. lA6l where we show the effects of appar- 
ent magnitude limits in the R and K-bands on the predicted redshift 
distributions for S'„(24^m) > 83piJy. In this plot, the redshift 
distributions are plotted as number per solid angle, without nor- 
malizing to unit area under the curve. The l eft and right panels re- 
spe ctively have th e redshift distributions of [Perez-Gonzalez et al. I 
and ICaputi et al. I overplotted. We concentrate on the comparison 
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Figure A2. Predicted galaxy redshift distributions in the four IRAC bands, for galaxies brighter than S v = 10/iJy. (a) 3.6 fim, (b) 4.5 (im, (c) 5.8 (im, and 
(d) 8.0 fim. The model curves (which all include the effects of dust) are as follows - blue: total; red: ongoing bursts; green: quiescent galaxies. The curves are 
normalized to unit area under the curve for the total counts. The median (250) and 10- and 90-percentile (zip , Z90) redshifts for the tota l counts in each band 
are also given in each panel. For 3.6 fim, the model predictions are compared with observational data from Franceschini et «T] 120061) (black dashed line), 
normalized to unit area as for the models. The error bars plotted on the observational data include Poisson errors only. 



with lCaputi et al. L since this has the simpler sample selection and 
more accurate redshifts. The model pred iction for K < 21.5 
(which is the magnitude limit used by ICaputi et al. I) is shown 
by the short-dashed blue line, while the prediction with no limit 
on the K-magnitude is shown by the solid blue line. The model 
d N/dz with no l imit on the K magnitude is most discrepant with 
the ICaputi et al. I data at z ~ 2, where it predicts ~ 2 times too 
many galaxies. This is directly related to the fact that the pre- 
dicted luminosity function at 2 = 2 at rest-frame wavelength 
8 /im (corresponding to observed wavelength 24 /im) and lumi- 
nosity ~ 10 1 1 Lp ) is also ~ 2 times too high compared to what 
ICaputi et al. I estimate from their data, as shown in Fig.[5jc). When 
the effect of the K < 21.5 limit is included, the predicted redshift 
distribution is closer to the observational data, but only 58% of the 
model galaxies are brighter than this K-band ma gnitude limit, as 
against 95% in the observed sample of ICaputi et al. I We conclude 
that the main reason for the discrepancy between the predicted and 
observed redshift distributions at 24 fim is that the model predicts 



a rest-frame 8 fim luminosity function at 2 ~ 2 which is somewhat 
too high at luminosities ~ 10 11 Lq, even though it reproduces quite 
well the general features of the evolution of the mid-IR luminosity 
function. 



24 Lacey et al. 



(a) 

\= 3.60 Aim 

S„= 1.0e-07 Jy 

S„= 1.0e-06 Jy 
S„= 1.0e-05 Jy 
S„= 1.0e-04 Jy 

S„= 1.0e-03 Jy 





4.50 Atm 

S„= 1.0e-07 Jy 

S„= 1.0e-06 Jy 
S„= 1.0e-05 Jy 
S„= 1.0e-04 Jy 

S„= 1.0e-03 Jy 




(c) 

A= 5.80 Aim 

S„= 1.0e-07 Jy 

S„= 1.0e-06 Jy 
S„= 1.0e-05 Jy 
S„= 1.0e-04 Jy 

S„= 1.0e-03 Jy 




1 1 1 1 1 I 1 1 r " 
(d) 

X= 8.00 Aim 

S„= 1.0e-07 Jy 

S„= 1.0e-06 Jy 

S„= 1.0e-05 Jy 

S„= 1.0e-04 Jy 

S„= 1.0e-03 Jy 



Figure A3. Predicted galaxy redshift distributions in the four IRAC bands, for different fluxes, (a) 3.6 /im, (b) 4.5 fim, (c) 5.8 Aim, and (d) 8.0 Aim. In this 
figure, the redshift distributions are for galaxies at a particular flux. Predictions are shown for fluxes S v = 0.1, 1, 10, 100 and 1000 fiJy, as shown in the key. 
In all cases, the model curves are normalized to unit area, and include the effects of dust. 
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Figure A4. Predicted galaxy redshift distributions in the three MIPS bands, 
for galaxies brighter than a specified flux, (a) 24 fim, S v > 83/iJy, (b) 70 
/im, S„ > lOmJy, and (c) 160 (im, S„ > lOOmJy. The model curves 
are as follows - blue: total; red: ongoing bursts; green: quiescent galax- 
ies. The curves are normalized to unit area under the curve for the total 
counts. The median (Z50) and 10- and 90-percentile (210, 290) redshifts for 
the total counts in each band are also given in each panel. F or 24 (im, the 
model p redictions are compar ed with observational data fromlCaputi et al. I 
1 2006a) (solid black line) and lPerez-Gonzalez et ai\ J2005l) (dashed black 
line), normalized to unit area as for the models. Th e error bars pl otted on 
the observational data include Poisson err ors only for | Caputi et al \ but also 
include errors in photometric redshifts for lPerez-GonzalezTTqZT 



Figure A5. Predicted galaxy redshift distributions in the three MIPS bands, 
for different fluxes, (a) 24 fim, (b) 70 jim, and (c) 160 fj,m. In this figure, 
the redshift distributions are for galaxies at a particular flux, as shown in the 
key in each panel. In all cases, the model curves are normalized to unit area, 
and include the effects of dust. 
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Figure A6. Predicted redshift distributions at 24fim, showing the effects of optical or near-IR magnitude limits. Model galaxies are selected with S„ > 83/^Jy 
together with the optical/NIR magnitude limits as shown i n the key. The fraction of 24 fi m sources brighter than each magnitud e limit is also g iven, (a) R-band 
magnitude limit. The observed redshift distribution from Perez-Gonzal ez et ai\ J2005t) is overplotted in black. No t elLe Floc'h et al. (2005) used R < 24 
and obtained 54% completeness, (b) K-band magnitude limit. The observed redshift distribution from lCaputi et al. 1 J2006al) (with K < 21.5) is overplotted. 
Magnitudes are on the Vega system. 



